home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / multiroots / dnewton.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.5 KB  |  187 lines
  1. /* multiroots/dnewton.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21.  
  22. #include <stddef.h>
  23. #include <stdlib.h>
  24. #include <stdio.h>
  25. #include <math.h>
  26. #include <float.h>
  27.  
  28. #include <gsl/gsl_math.h>
  29. #include <gsl/gsl_errno.h>
  30. #include <gsl/gsl_multiroots.h>
  31. #include <gsl/gsl_linalg.h>
  32.  
  33. /* Newton method using a finite difference approximation to the jacobian.
  34.    The derivatives are estimated using a step size of 
  35.  
  36.    h_i = sqrt(DBL_EPSILON) * x_i 
  37.  
  38.    */
  39.  
  40. typedef struct
  41.   {
  42.     gsl_matrix * J;
  43.     gsl_matrix * lu;
  44.     gsl_permutation * permutation;
  45.   }
  46. dnewton_state_t;
  47.  
  48. static int dnewton_alloc (void * vstate, size_t n);
  49. static int dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
  50. static int dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx);
  51. static void dnewton_free (void * vstate);
  52.  
  53. static int
  54. dnewton_alloc (void * vstate, size_t n)
  55. {
  56.   dnewton_state_t * state = (dnewton_state_t *) vstate;
  57.   gsl_permutation * p;
  58.   gsl_matrix * m, * J;
  59.  
  60.   m = gsl_matrix_calloc (n,n);
  61.   
  62.   if (m == 0) 
  63.     {
  64.       GSL_ERROR_VAL ("failed to allocate space for lu", GSL_ENOMEM, 0);
  65.     }
  66.  
  67.   state->lu = m ;
  68.  
  69.   p = gsl_permutation_calloc (n);
  70.  
  71.   if (p == 0)
  72.     {
  73.       gsl_matrix_free(m);
  74.  
  75.       GSL_ERROR_VAL ("failed to allocate space for permutation", GSL_ENOMEM, 0);
  76.     }
  77.  
  78.   state->permutation = p ;
  79.  
  80.   J = gsl_matrix_calloc (n,n);
  81.  
  82.   if (J == 0)
  83.     {
  84.       gsl_permutation_free(p);
  85.       gsl_matrix_free(m);
  86.  
  87.       GSL_ERROR_VAL ("failed to allocate space for d", GSL_ENOMEM, 0);
  88.     }
  89.  
  90.   state->J = J;
  91.  
  92.   return GSL_SUCCESS;
  93. }
  94.  
  95. static int
  96. dnewton_set (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
  97. {
  98.   dnewton_state_t * state = (dnewton_state_t *) vstate;
  99.   size_t i, n = function->n ;
  100.   int status;
  101.  
  102.   status = GSL_MULTIROOT_FN_EVAL (function, x, f);
  103.   if (status)
  104.     return status;
  105.  
  106.   status = gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON,
  107.       state->J);
  108.   if (status)
  109.     return status;
  110.  
  111.   for (i = 0; i < n; i++)
  112.     {
  113.       gsl_vector_set (dx, i, 0.0);
  114.     }
  115.  
  116.   return GSL_SUCCESS;
  117. }
  118.  
  119. static int
  120. dnewton_iterate (void * vstate, gsl_multiroot_function * function, gsl_vector * x, gsl_vector * f, gsl_vector * dx)
  121. {
  122.   dnewton_state_t * state = (dnewton_state_t *) vstate;
  123.   
  124.   int signum ;
  125.  
  126.   size_t i;
  127.  
  128.   size_t n = function->n ;
  129.  
  130.   gsl_matrix_memcpy (state->lu, state->J);
  131.  
  132.   {
  133.     int status = gsl_linalg_LU_decomp (state->lu, state->permutation, &signum);
  134.     if (status)
  135.       return status;
  136.   }
  137.   
  138.   {
  139.     int status = gsl_linalg_LU_solve (state->lu, state->permutation, f, dx);
  140.     if (status)
  141.       return status;
  142.   }
  143.  
  144.   for (i = 0; i < n; i++)
  145.     {
  146.       double e = gsl_vector_get (dx, i);
  147.       double y = gsl_vector_get (x, i);
  148.       gsl_vector_set (dx, i, -e);
  149.       gsl_vector_set (x, i, y - e);
  150.     }
  151.   
  152.   {
  153.     int status = GSL_MULTIROOT_FN_EVAL (function, x, f);
  154.  
  155.     if (status != GSL_SUCCESS) 
  156.       {
  157.         return GSL_EBADFUNC;
  158.       }
  159.   }
  160.  
  161.   gsl_multiroot_fdjacobian (function, x, f, GSL_SQRT_DBL_EPSILON, state->J);
  162.  
  163.   return GSL_SUCCESS;
  164. }
  165.  
  166.  
  167. static void
  168. dnewton_free (void * vstate)
  169. {
  170.   dnewton_state_t * state = (dnewton_state_t *) vstate;
  171.  
  172.   gsl_matrix_free(state->J);
  173.   gsl_matrix_free(state->lu);
  174.   gsl_permutation_free(state->permutation);
  175. }
  176.  
  177.  
  178. static const gsl_multiroot_fsolver_type dnewton_type =
  179. {"dnewton",                /* name */
  180.  sizeof (dnewton_state_t),
  181.  &dnewton_alloc,
  182.  &dnewton_set,
  183.  &dnewton_iterate,
  184.  &dnewton_free};
  185.  
  186. const gsl_multiroot_fsolver_type  * gsl_multiroot_fsolver_dnewton = &dnewton_type;
  187.